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Abstract. We study the dynamics of a one-dimensional Bose gas after a sudden change of 
the interaction strength from zero to a finite value using the numerical time-evolving block 
decimation (TEBD) algorithm. It is shown that despite the integrability of the system, local 
quantities such as the two-particle correlation g( 2 \x,x) attain steady state values in a short 
characteristic time inversely proportional to the Tonks parameter 7 and the square of the 
density. The asymptotic values are very close to those of a finite temperature grand canonical 
ensemble with a local temperature corresponding to initial energy and density. Non-local 
density-density correlations on the other hand approach a steady state on a much larger time 
scale determined by the finite propagation velocity of oscillatory correlation waves. 
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1. Introduction 

The dynamics of interacting quantum systems from an initial non-equilibrium state constitutes 
a major challenge for many-body theory. In particular the question of thermalisation of 
integrable models regained attention recently due to the experimental progress in ultra-cold 
gases. As demonstrated in a beautiful experiment by Kinoshita et al. HI for the example 
of a ID ultra-cold Bose gas, integrable systems do not fhermalise in the usual sense, i.e., 
reduced density-matrices relax on considerably different time scales than in the absence of 
integrability. 

The speciality of the relaxation dynamics of integrable systems has been attributed to the 
presence of an infinite set of constants of motion with local character, i.e., which can be written 
as sums of operators acting only over a finite spatial range. Although thermalisation has been 
studied in a large body of theoretical papers it remains a largely unsolved problem. Most 
studies of specific models have been done either for non-interacting particles or systems 
that can directly be mapped to free systems such as hard-core bosons [3], the Luttinger model 

01, or the 1/r fermionic Hubbard model 0. 

In the present letter we analyse the dynamics of a ID Bose gas with s-wave scattering 
interactions, described by the Lieb-Liniger (LL) model, after a sudden quench of the 
interaction strength from zero to a finite value, covering the full range from weak to strong 
interactions. Performing numerical simulations using the time evolving block decimation 
algorithm (TEBD) (6113, we show that local quantities, in particular the local two-particle 
correlation g^ 2 \0, 0; t), do approach steady-state values on a short time scale determined only 
by the Tonks parameter 7 and the particle density q. This shows that although non-local 
quantities such as the momentum distribution do not approach a steady state over long times 
jH, there is an equilibration in a local sense. Furthermore the asymptotic values of g^ 2 \x, x) 
are very close to those obtained from a thermal Gibbs ensemble 10, with temperature and 
chemical potential determined by the initial conditions and the amplitude of the interaction 
quench. Thus it is possible to define local temperature and chemical potential and the 
influence of constants of motion other than total energy and particle number is very small, if 
present at all. Non-local quantities such as the density-density correlation approach a steady- 
state distribution on a larger time scale by way of correlation waves propagating out of the 
sample. 

2. LL model and lattice approximation 

A Bose gas in one spatial dimension is described by the Hamiltonian 



in units were h — m — k-Q — 1. Here is the field operator of the Bose gas in second 
quantisation, V (x) some possible trap potential, and g the strength of the local particle-particle 
interaction. The latter is characterised by the dimensionless Tonks parameter 7 = g/g, where 




^(2) + 
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q = (W(x)^f(x)) is the ID density of particles. Specifically we consider here a system 



initially prepared in the non-interacting ground state. 

The initial canonical state has locally only diagonal elements. In the course of 
interactions non-diagonal elements are not created. Thus the reduced local density matrix 
is entirely determined by the number distribution, and the quantities of interest are the density 
p and the local two-particle correlation g^ (x, x, t), where 



In principle also the corresponding higher-order moments are nonzero. They will, however, 
not be considered here. 

We study the dynamics after the interaction quench numerically by means of the the 
time-evolving block decimation algorithm J6l [Vj . As for any numerical methods this requires 
the use of a discretised version of the Lib-Liniger model. The discretization of the LL model 
leads to the (non-integrable) Bose-Hubbard model [[Toll 



Here J = 1/ (2Ax 2 ), U = g/Ax, and D = 1/ Ax 2 , with Ax being the lattice constant of the 
discretization grid. The appropriateness of discretised lattice models to describe continuous 
interacting Bose or Fermi gases in the limit Ax — > has been discussed and verified in a 
number of earlier papers ifTOlfTTTl . Note furthermore that using the boson-fermion duality in 
ID the LL model can be mapped to an integrable lattice model, the spin 1/2 XXZ model ifTTTl . 
For some data sets we used both, Equation ©, and the XXZ lattice model to verify that in the 
considered limit the non-integrability of © has no influence on the results. 

It turned out that for the simulation of dynamics the necessary grid sizes are much smaller 
than for equilibrium simulations ifTOl [HJ. Empirically we found that in order to minimise 
lattice artifacts resulting into numerical errors, the average number of particles per lattice site 
(n) should be small compared to where j — g/g corresponds to the density at the centre 
of the cloud, i.e. g — > g(x = 0). This can be explained as follows: The interaction energy 
of a two-particle collision in the lattice, i.e., gj Ax should be smaller than the bandwidth 
of the lowest Bloch band, ~ I/Ax 2 , in order not to see lattice artifacts. At the edges of the 
cloud where collisions become less probable the condition is less strong. To accommodate the 
requirement of a very small (n) at the centre of the cloud we used space dependent grid sizes 
such that the average boson number per site was constant for the centre part of the particle 
distribution. Nevertheless to approximate the continuous model sufficiently well, very fine 
grids were needed leading to rather large lattice sizes on the order of up to 2880, which is 
rather challenging. In order to illustrate the effects of discretization we have plotted in Figure 
[1^ the local two-particle correlation (see following section) for 7 = 200/9 and increasing 
lattice sizes L, corresponding to finer grids. One clearly recognises oscillation artifacts which 
only slowly disappear with increasing L. 

The convergence of the TEBD scheme was checked by varying the bond dimension 
X of the matrix product state (MPS) and calculating the truncation error in the state norm 
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g{x)g{y) 



(2) 




Xj = j Ax. (3) 
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Figure 1. a) Time evolution of normalised local two-particle correlation g( 2 )(0, 0;i) (see 
section^ for N = 9 particles and 7 = 200/9 for increasing lattice length, corresponding to 
finer grid sizes. One clearly recognises oscillations which are lattice artifacts and which only 
disappear for the largest lattice sizes, b) Accumulated truncation error of the norm of the MPS 
in the dynamical TEBD algorithm for 7 = 200/9 and increasing bond dimension \- 

accumulated during the time evolution. In Figure [Tt> the accumulated truncation error is 
plotted for 7 = 200/9 and increasing values of x from 25 to 200. One recognises that for 
the maximum value of x — 200 which we were able to use, the truncation error is below the 
level of 10~ 3 for the time scale of interest. This value is larger then the accuracy typically 
reached in ground state calculations. However we are not at the point where the cut-off 
explodes, which typically happens in dynamical simulations at some point. Finally the matrix 
dimension required to achieve a given accuracy does not depend on the discretization length, 
i.e. the number of lattice sites used. It is rather the number of particles which determines the 
complexity of the calculations. Thus the restriction to a moderate particle number allows us 
to work on lattice large compared to other applications of the algorithm. 
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3. Local Dynamics 

In order to be able to perform numerical simulations with a fixed number of particles (up to 1 8, 
which corresponds to the experiments in references lfT2l[T3l ) we have to work with a finite size 
system. Therefore we assumed an initial weak harmonic trapping potential V(x) = \u 2 x 2 . 
Open or periodic boundary conditions, which can also be dealt with by the TEBD algorithm, 
could be used alternatively. Initially the Bose gas is in the canonical ground state (T = 0) of 
non-interacting bosons in the trap, for which the matrix product representation is analytically 
known, since it is a product of single particle states. At t = we suddenly switch the 
interaction strength from zero to a finite value. At the same time the trap, the only purpose 
of which is was the preparation of an appropriate initial state, is switched off. On the time 
scales we are interested in, the density distribution does not change, so that the presence of a 
trap would be of no relevance. This also allows to apply the results of the present analysis to 
a homogeneous gas. 

The initial state has a Gaussian density distribution q{x) = N ^T^e'^ with / OS( - 



being the oscillator length. Figure [2] shows the time evolution of g( 2 \0,0;t) with time 
normalised to a characteristic time scale t ia for different values of 7. One recognises after 
an initial phase a power-law decay with an exponent that is monotonous in the interaction 
parameter. At times close to t- m a steady state value is attained indicating that a local 
equilibrium is reached. I.e., although globally a LL gas does not thermalise HI, local 
quantities do. The time scale t i£L of the local dynamics can be estimated from Equation ©. The 
repulsive interaction Un(n — 1) causes particle number fluctuations to be driven out of a given 
lattice site. This happens in the following way: Initially all components of the state vector 
have the same phase and tunnelling has no effect. However, due to the interaction, components 
with different particle number attain a differential phase shift and are subsequently coupled to 
states in adjacent lattice sites by tunnelling with rate J. Since in the limit Ax — > we have 
J 3> U, the maximum rate of this process is limited by the average interaction energy per 
particle U(n). Thus we have 

t- = — = - = — (4) 

^la ttI"\ 2' VV 

U (n) gg 7^ 

Note that already for moderate interaction strength, 7 ^> 1/N 2 , this time is much shorter as 
for example the oscillation time t osc in the trap. 

/2 sc 1 

Note furthermore that although the characteristic time of the expansion of the gas after 
switching off the trap becomes much shorter than t osc for larger interactions, we found that 
the density profile did not change on the scale of t[ a even for the largest values of 7 used, since 
also t ia ~ l/7.Note furthermore that although the characteristic time of the expansion of the 
gas after switching off the trap becomes much shorter than t osc for larger interactions, it will 
be large compared to t ia . This is because the kinetic energy transferred to the particles will be 
of the order 7 and therefor their characteristic speed will be of the order ^7 only. Accordingly 
we found numerically that the density profile did not change on the timescale £ ia even for the 



Fermionisation dynamics of a strongly interacting ID Bose gas 



6 




Figure 2. (Colour online) Time evolution of normalised local two-particle correlation 
</ 2 ) (0, 0; t) after a sudden switch on of interactions at t = obtained from a numerical TEBD 
simulation for 9 particles initially prepared in the non-interacting ground state of a harmonic 
trap. An intermediate power-law decay with an exponent that is monotonous in 7 is apparent. 
The lattice size was up to L = 2880 for the strongest interactions corresponding to a lattice 
spacing of Ax/l osc w 6.15 • 1CP 4 at the trap centre. 



largest values of 7 used. Whether or not the thermalised local correlation will adiabatically 
follow the density evolution after longer times, i.e. when the expansion of the cloud sets in, 
cannot be concluded from our simulations. We would however expect such a behaviour. 

The fluctuations in the plots are artifacts of the discretization, which leads to an 
oscillatory behaviour of on top of the continuous-system time evolution. These artifacts, 
which are most pronounced for larger interactions, could not be eliminated completely even 
for the smallest grid sizes used. As a result the asymptotic values of <7^(0, 0, t) can only be 
given with a certain error. 

In figure [3] we have plotted the exponents obtained from a fit to the curves in Figure [2] 
which for intermediate times follows a power law 

0< 2 >(O,O;f)~ (IX . (6) 

The exponent is a monotonous function of the interaction strength and slowly approaches the 
limit — 1 for 7 -> 00, i.e., for a Tonks-Girardeau gas lfl4l . 

We now want to analyse the local state of the system in the stationary limit. In particular 
we will show that the local steady-state can be well described by the usual finite-temperature 

(2) 

Gibbs state. To this end we calculate the expected asymptotic value <7yy(0, 0) fr° m me 
thermodynamic Bethe Ansatz 0. The system is initially prepared in its non-interacting 
ground state, so we have <7 ini (0, 0) = 1 — 1/N, which in the thermodynamic limit — > 00 
approaches unity. The energy of this state with respect to the non-interacting Hamiltonian is 
0. At time t — the interaction is switched to a finite strength g > and the expectation value 
of the interaction energy immediately after the quench is given by 

£ ia = [dx^-(¥\x)^ 2 (x))= fdx?-g ( - 2 \x,x)Q 3 (x). (7) 
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Figure 3. (Colour online) Exponents a — 1 of the intermediate power-law decay of g 2 (0, 0; t) 
in figure[2]as function of 7. Error-bars indicate systematic fitting error. 



Since in a homogeneous system there is no x-dependence the energy per particle is 



E 

N 



t=o+ 



N 



1%. 



(8) 



t=o+ 



Here we have introduced the critical temperature T c in one dimension T c = g 2 /2. One 
recognises that the system is in a highly excited non-equilibrium state after the quench if 
7 > 1. Using the energy per particle, the density g and the Tonks parameter 7 as input 
parameter, we can extract a temperature T of a corresponding thermal Gibbs state by inverting 
the Yang- Yang equations of the thermodynamic Bethe Ansatz [|9]| for the excitation-energy 
and particle densities e(q),n(q) in momentum space 

\ 2 rp poo 

6(A) = y - ^ - ^ / 5 K(X t In (1 + e-«V/T) , 



Here K(\,£) 



29 



-(A-€) : 



and q = J dAn(A). With the help of the Hellmann-Feynman 



theorem we can then obtain the value ^yy(O'O) corresponding to the Gibbs state at 
temperature T lfT5l : 



^ (0 '° ) = ?^ /(7 ' T) 
with 7(7, T) being the free energy per particle 

T 



2 d 



/(7,r) = /* 



2vrp 



/oo 
ln(l + e-^/ T ) 
-oo 



In the limit 1 < r < 7 2 , where r = T/T c , © attains the simple form 031 

2r 
7' 



S? Y (0,0) 



(9) 



(10) 



(11) 
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Figure 4. (Colour online) Steady-state values of (0, 0) (left scale) for different values 
of the interaction parameter 7 after the interaction quench obtained from TEBD simulations 
using 9 (red circles, corresponding to figure [2]i and 18 (green triangles) bosons. Solid line: 
value from thermal Gibbs state in the thermodynamic limit; dot dashed line: temperature (right 
scale) corresponding to the given energy per particle in the thermodynamic limit. The error 
bar reflects discretization error estimated by comparing steady-state values for L = 720 and 
L = 1440 lattice sites as well as error resulting from finite MPS matrix dimension obtained 
from comparing results for bond dimension \ = 100 and 200. 

In Figure |4] we have plotted the values of <7yy(0, 0) from the thermal Gibbs state in the 
thermodynamic limit as function of the interaction strength 7 (solid line). Also shown are the 
steady-state values obtained from the numerical simulation in Figure|2] The error bar indicates 
uncertainties which are here due to discretization artifacts and error estimates obtained from 
comparing simulations with MPS bond dimensions \ — 100 an d 200. It is available only for 
one parameter set, since the variation of the discretization length and the bond dimension is 
numerically expensive. However we expect it to be of about the same relative size for all data 
points. One recognises that g^ 2 \0, 0; t) attains in the long-time limit values which are close 
to that of the thermal Gibbs state. One should note that the steady-state values for the largest 
values of 7 (7 = 500/9 and 18 bosons as well as 7 = 1000/9 and 9 bosons) are slightly 
overestimated in the simulation due to the remaining grid artifacts since here (77)7 ~ 0.73 is 
no longer small compared to unity. Also shown is the asymptotic local temperature of the gas 
in units of the degeneracy temperature. For large values of 7, T 3> T c , i.e., after relaxation 
the gas is in a state with large local temperature. 

4. Non-local Relaxation 

We now discuss the dynamics of non-local quantities. Specifically we consider the non-local 
two-particle correlation g^(Q,x;t). In Figure |5] we have plotted g( 2 \0,x;t) for different 
times after the interaction quench. One recognises that while the local correlations attain 
a steady-state value on a short time scale, the non-local evolution happens much slower. 
Switching on the particle-particle repulsion leads to a fast reduction of the probability to find 
two particles at the same position. Associated with this is a correlation flow to larger distances 
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Figure 5. (Colour online) Time evolution of non-local density-density correlations 
<7^(0, x; t) for 7 = 200/9. x = denotes the centre of the cloud. One recognises the 
formation of expanding correlation waves. The dashed blue line shows the approximation [\3[ 
to the non-local correlation in a thermal Gibbs state from [16) multiplied by </ 2 )(0,0,i = 
0) = 8/9 to account for the final particle number (N = 9) used in the simulation. 



leading to expanding correlation waves. For very short times the propagation velocity of 
correlation waves is faster than the Fermi velocity vf = ttq- But at the largest time shown in 
Figure [5] corresponding to t = 0.01t osc , the maximum of the correlation wave has travelled a 
distance of approximately Ax = 0.12/ osc which is consistent with the speed of sound which 
for large values of 7 approaches 

v s = v F (l-^\. (12) 

The buildup of a maximum that behaves like a wavefront can be understood as follows: 
The integral over space of 5^(0, x; t) is (in a homogeneous system) a constant with respect 
to time due to particle number conservation. The numerator of g^(0, x; t) is proportional to 
the joint probability distribution to find a particle at position x given that there is one particle 
at the origin. 

So as the quench can not change (0, x; t) significantly outside the light cone given by 
the Fermi velocity, the reduction of the probability to find two particles close together must 
be accompanied by an increase at finite distance. 

As one can see from Figure \5\ also the non-local correlation function approaches at least 
for smaller distances in the large-time limit that of the thermal Gibbs state with temperature 
and density given by the initial conditions and the Tonks parameter 7. For comparison we 
have plotted an approximation to the finite-temperature non-local g( 2 ' from reference lfT6l 
which holds in the regime 1 < r < 7 2 



g®(0,x) = l- 




TIT I X 

7 v> 



e -M*M\ (13) 



Here A r = ^/Att/tq 2 is the thermal de Broglie length. 
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5. Experimental Observation 

In the following we discuss the possibility to test the local relaxation in an experiment. For 
this we make use of the fact that by energy conservation the interaction energy lost by the 
decrease of g^(Q, 0) must be gained as kinetic energy, E khl (t) = E- 1& (t = 0) — E ia (t) and the 
kinetic energy therefore directly gives the value of «7^(0, 0, t) in the homogeneous case: 

£ kin (t) = J dx | (1 - s (2) (0,0; t )) (?. (14) 

If the interaction is turned on at t = and turned off abruptly at some time t = t\ the kinetic 
energy is the only remaining in the system and can be used to measure g^ 2 \Q, 0; t\). 

o 771 final 

^)(0,0;t 1 ) = l-4^r- (15) 

7£T JN 

In an the experimental setup, the gas must be confined e.g. by an harmonic trapping potential. 
So the initial non-interacting state has a Gaussian density distribution. It is also a good 
assumption, that the correlations decay locally as in the homogeneous system corresponding 
to the local density provided the density q{x) remains constant over the time scale of interest. 
This is indeed the case, if t- mt = l/(7p 2 ) <C l 0&c / v s ~ fcsc/p This means, that the Tonks 
parameter must be large compared to -k, which is of course the case we are interested in. We 
note that the region in the wings of the density distribution which does not fulfil this constraint 
gives a negligible contribution to the total interaction energy. Of course measuring the kinetic 
energy in the trap gives only an average of 9 Jf^ over the trap. 

6. Summary 

Using the time-evolving block decimation scheme we have numerically analysed the 
dynamics of a ID Bose gas (LL-model), after an interaction quench from zero to a finite value. 
Although globally the ID Bose gas does not fhermalise, we have shown that local quantities 
attain a steady-state value on a time scale £ ia = ( , yg 2 )~ 1 . Within the achievable accuracy these 
values are consistent with the assumption that local quantities relax to a thermal Gibbs state 
with local temperature determined by the initial energy and chemical potential. Non-local 
quantities such as the density-density correlation relax on a much longer time scale set by the 
velocity of sound by means of correlation waves propagating out of the sample. 
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